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SUMMARY 


The transport of a reacting permeant diffusing through a thin membrane is extended to more realistic 
dissociation models. A new nonlinear analysis of the reaction-diffusion equations, using implicit finite- 
difference methods and direct block solvers, is used to study the limits of linearized and equilibrium 
theories. Computed curves of molecular oxygen permeating through hemoglobin solution are used to 
illustrate higher-order reaction models, the effect of concentration boundary layers at the membrane 
interfaces, and the transient buildup of oxygen flux. 


INTRODUCTION 


Membranes have been found to be very useful for separating both gaseous and aqueous mixtures. 
Artificial membranes with very selective separation properties are now commercially available and may 
find many applications in advanced life support systems for future space missions. The concentration 
and separation of aqueous wastes have been studied extensively using the reverse-osmosis process, 
which depends on the passive action of membrane systems (Lonsdale, 1986). Another, and very intrigu- 
ing, aspect of membrane separation is the facilitated transport mechanism (Schultz, 1986). This mode of 
diffusive transport was discovered in experimental studies of oxygen diffusion through films of hemo- 
globin solution, and has been the subject of many studies. Facilitated transport depends on the action of 
a reversible chemical reaction to increase the flux of a permeant many times over that existing in a pas- 
sive medium. This is a particularly effective method for the selective transport of specific molecules 
such as oxygen through a liquid membrane of hemoglobin in blood plasma. This paper presents a new 
analysis of facilitated transport in a particular biological process — the enhanced diffusion of oxygen 
using oxygenated hemoglobin as a carrier. This analysis should also be useful in the development of new 
carrier species, both natural and synthetic, for the efficient separation and concentration of gaseous com- 
ponents of life support systems (such as carbon dioxide) for long-duration space missions. 

Hemmingsen and Scholander (1960) and Scholander (1960) demonstrated that the oxygen flux 
through thin liquid films is greatly enhanced by the presence of the protein hemoglobin. Furthermore, 
the degree of enhancement was shown to be sensitive to the oxygen concentration difference across the 
membrane. In an attempt to explain these important observations, many mathematical models were 
developed (Collins, 1960; Wang, 1961; Fatt and La Force, 1965; Friedlander, 1965). These early papers 
made the critical assumption that oxygen and hemoglobin were in chemical equilibrium across the film. 
This is a reasonable assumption for the oxygen-hemoglobin system for films thicker than about 
100 microns. In the equilibrium state, the degree of facilitation can be simply computed with algebraic 
techniques or even estimated qualitatively using the oxygen dissociation curve. 

The steady-state diffusion of oxygen through hemoglobin solution was computed for the case of 
nonequilibrium chemical reactions by Kutchai et al. (1970). The equations governing the process are a 
system of coupled nonlinear ordinary differential equations. Using a numerical algorithm based on the 
quasi-linearization scheme, they were able to predict the overall facilitation as well as the detailed distri- 
bution of oxygen and hemoglobin across the liquid membrane. The system of equations displays a 
boundary layer effect at the air/film interfaces which has been the source of past difficulties in the 



analysis of these systems. In Kutchai et al. (1970), a special iteration scheme was developed to obtain 
converged solutions for increasing film thicknesses. An important conclusion of this work was that in 
vitro studies of facilitated transport are adequately predicted by equilibrium theories, but for systems of 
the dimensions of human erythrocytes, the more complex nonequilibrium analysis must be used since 
the effective film thickness cannot be accommodated by equilibrium theory. 

Following these studies, a number of new analysis methods were introduced. A numerical scheme 
based on the concept of “local linearization” was introduced by Gonzalez-Femandez and Atta (1981, 
1982). An approach based on the finite element technique was used by Jain and Schultz (1982). This 
latter work, in particular, includes a very comprehensive discussion of the boundary-layer effect and 
introduces one particular method to treat this difficulty. On the analytical side, this problem has been 
treated as an example of a “singular perturbation” (Goddard et al., 1970; Schultz et al., 1974) where the 
governing differential equations are such that the highest derivative is multiplied by a small parameter. 

In the limiting case of an infinitely thick membrane, this parameter vanishes and the system is in chemi- 
cal equilibrium. In this limit, it is not possible to satisfy all the boundary conditions, and this has been 
the source of some confusion in the past. This boundary-layer effect is also the basis for the term “stiff- 
ness” that has been used for this class of transport processes and it causes many difficulties in the solu- 
tion of facilitated transport problems. (It is interesting that this is a relatively benign singular pertur- 
bation since the derivative of the unknown function is given as a boundary condition. Its limiting value, 
however, can easily be computed from equilibrium theory. In the original singular perturbation problem 
concerning the flow of a viscous fluid, the value of the function is given and its derivative is unknown. 
The limiting case is the nonphysical slip-flow solution with the derivative undefined.) All of the work 
mentioned above treats the steady-state transport and assumes that the chemical reaction is a second- 
order process. 

The flow of oxygen across a thin film is not an instantaneous process and, in fact, the development 
time is a strong function of the system parameters. The study of unsteady facilitated transport is not as 
well developed as the steady problem. Recent studies of unsteady facilitated transport by Ruckenstein 
(1982) and Varanasi (1983) also used the local linearization method mentioned above and assumed 
second-order chemical reactions. 

In this paper, the unsteady transport of oxygen across a hemoglobin solution with higher-order 
chemical reactions is studied using a finite-difference algorithm based on the method of time lineariza- 
tion. Before applying the algorithm, the governing systems of equations are recast into a first-order 
matrix-differential equation. This has the distinct advantage that all fluxes are treated as dependent vari- 
ables and are computed directly as the solution evolves. Another advantage is that the transport of mul- 
tiple permeant- species can be easily accommodated — the only change being the order of the matrices 
involved. A third advantage of the matrix approach is that the eigenvalue spectrum can be easily com- 
puted, which gives a good physical understanding of difficulties found in previous solution methods. It 
is shown that this problem is yet another example of a class of stiff problems that arise frequently in 
practical situations. The facilitated transport problem is stiff in the spatial sense, which means that 
eigenmodes of disparate scales are generated. 

It has been demonstrated in the theory of nonequilibrium high-speed gas dynamics (Lomax and 
Bailey, 1967) that implicit numerical methods are almost mandatory when solving stiff problems. In this 
paper, an implicit algorithm based on a noniterative block-triadiagonal solver is implemented in the 


2 



development of an efficient solution process. A linearized form of the reaction-diffusion equation is 
derived that highlights the disparate time scales. This linearized system is used to develop an explicit 
formula for the size of the boundary layer region. In the following sections, solutions to the full nonlin- 
ear systems are presented for the oxygen-hemoglobin system and comparisons made with previously 
published results. The effect of initial conditions and physically plausible oxygen dissociation curves on 
the final steady-state flux is computed using the system parameters of Scholander (1960), Hemmingsen 
and Scholander (1960), and Kutchai et al. (1970). 


THE OXYGEN DISSOCIATION CURVE AND FACILITATED TRANSPORT 


A very important measure of the oxygen affinity of hemoglobin is its characteristic dissociation 
curve. In most theoretical treatments, a second-order reaction is assumed which restricts the curve to one 
branch of a hyperbola. This is a good approximation for oxygen transport through myoglobin, but is not 
as accurate for the cooperative binding that occurs in hemoglobin molecules. A simple analytical 
approximation to the oxygen dissociation curve is (Stryer, 1981, p. 67): 



where Y is the ratio of oxyhemoglobin to the total hemoglobin present and a is the solubility coefficient 
(assuming Henry’s Law). p02 is the partial pressure of oxygen ( 0 tpO 2 is its concentration in moles/cc) 
and pO2,0.5 is the oxygen partial pressure when Y = 0.5, which is a commonly used reference point for 
the dissociation curve. The effect of the exponent n will be investigated for the cases n = 1 (hyperbolic 
dissociation curve) and n = 2.8 (Hill’s exponent, sigmoidal dissociation curve). 

Dissociation curves are shown in figure 1 for (1) an approximate second-order curve that matches 
the 50% saturation point in Scholander (1960), (2) an approximate 3.8th-order reaction that matches the 
same 50% saturation point, and (3) an approximate 3.8th-order reaction that matches the 3.8th-order 
reaction for hemoglobin at physiological conditions. The experimental data show a very rapid rise to 
100% saturation that is better matched to the sigmoidal shaped curve than the hyperbola. The increased 
oxygen affinity for the in vitro experiments referenced in Scholander (1960) may be due to the absence 
of certain substances in the intact erythrocytes. 

These curves can also be used to assess the relative degree of facilitation assuming equilibrium 
conditions throughout the film. The slope of a line from the origin to any point on the dissociation curve 
is proportional to the facilitation factor for the transport of oxygen from the indicated value of pC>2 to a 
vacuum. It is clear from figure 1 that for pC>2 < pC>2,o.5 the facilitation factor for a second-order reaction 
is increased relative to that using Hill’s exponent, while the reverse is true for pC >2 > p02,o.5- Similar 
comparisons can be made for other oxygen concentration differences. Further details of the facilitation 
process for nonequilibrium conditions must be made with the aid of the governing differential equations. 
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GOVERNING EQUATIONS AND LINEARIZED FORMULATION 


The continuity equation for the transient diffusion of a dilute permeant (oxygen) in a solvent that 
supports a reversible chemical reaction (hemoglobin solution) is generally expressed as a reaction- 
diffusion equation: 


^ = DV 2 C-R (2) 

dt 

where C is the species concentration, D is the diffusivity, and R is the rate of reaction of the species C. 
In the case of a one-reaction/three- species system in a one-dimensional membrane, equation 2 is 
replaced by three coupled partial differential equations: 

d[0 2 ] _ £a 3 2 [Q 2 ] 

3t ~ L 2 dx 2 R 

3[Hb] _ P B a 2 [Hb] (3) 

dt L 2 dx 2 

3[Hb0 2 ] D b a 2 [Hb0 2 ] 

— — — . — + R 

5t L 2 dx 2 

The bracketed quantities denote concentrations of oxygen, hemoglobin, and oxygenated hemoglobin, 
respectively. The diffusivity for hemoglobin and oxygenated hemoglobin are assumed equal. The main 
justification for this assumption is the much greater molecular weight of the hemoglobin protein com- 
pared with the oxygen molecule (Kutchai et al., 1970). The distance x across the membrane is normal- 
ized by the thickness L. The solution to equation 3 must satisfy the following boundary conditions at 
every instant of time: 


9[Hb0 2 ] 
[° 2 1 = Co fa 


3[Hb] 

dx 


@ x — — 0. 5 


[0 2 ] = C! 


3[Hb0 2 ] 

dx 


0[Hb] 

dx 


@ x = +0.5 


(3a) 


which state that under the influence of the driving concentration gradient (Co ~ Ci), the hemoglobin- 
containing species are trapped in the membrane. The most useful result of the analysis is the flux of the 
permeant, -D / \/L3[02]/3 x, at the membrane exit. 

Much of the literature on facilitated transport is concerned with the steady-state solution applied to 
the following idealized second-order chemical reaction: 
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0 2 + Hb <-> HbQ 2 


(4) 


This reaction will be compared with the reaction 

n0 2 + Hb HbO^ (5) 

which reproduces the experimental dissociation curve more accurately when n = 2.8 as shown as curve 2 
in figure 1. It should be noted that the factor n is just an artifact used to match experimental data and is 
not an actual stoichiometric coefficient. 

The kinetic constants ki and k2 governing the rate of reaction are 

R = kj [0 2 f[Hb] - k 2 [Hb0 2n ] (6) 

and the value of n is taken as either 1 or 2.8. The rate of reaction is a nonlinear function of the species 
concentration, and this is the main source of difficulty in solving the system of equations analytically. 

With the assumption that the diffusivities of the hemoglobin-containing species are equal, which is 
quite good for this system, the relation between the Hb-containing molecules is simply (Kutchai et al., 
1970) 


[Hb] + [HbOj = Q* (7) 

where Q* is the total hemoglobin concentration (a constant in moles/cc). 

Aside from the nonlinearity, the other major difficulty is the disparate time scale in equation 3 
between the diffusive process (the space derivative term) and the chemical reaction (the forcing 
function). This important effect can best be illustrated with the aid of a linear model and an appropriate 
scaling. A common assumption used in the literature to linearize the facilitated transport problem (Noble 
et al., 1986) is that [Hb] is always at a constant level in equilibrium with [O 2 ] and[Hb02]. In the case of 
n=l: 


[Hb] = 


k 2 [Hb0 2 ] 

kj[0 2 ] 


( 8 ) 


Eliminating the hemoglobin-oxygen complex from equations 7 and 8 yields the final formula for 
[Hb]: 


[Hb] = 


Q* 

1 + K[0 2 ] 


(9) 
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where K is the equilibrium constant ki/k2- The oxygen concentration is some representative value such 
as the mean boundary concentration Ao- With one of the reacting species constant, the system of equa- 
tions reduces to a coupled set of linear ordinary differential equations: 


Da 0 
0 D b 


a 2 o 2 


kiQ* 

3 x 2 


1 + KAo 

a 2 [Hb 0 2 ] 


-kiQ* 

ax 2 


1 + kA 0 



0 


0 


( 10 ) 


For further analysis, it is useful to consider the system as a set of first-order ordinary differential 
equations. Let the vector^ be defined as the nondimensional quantities {[O2VQ*, [C> 2 ] 7 Q*, [Hb02]/Q*, 
[HbC>2]YQ* }T. The final 4x4 system of linear equations with constant coefficients is: 


w i 


0 

1 0 0" 


Wi 

w 2 
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*3 

Ti 

0 tt 0 
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w 2 

w 3 


0 
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w 3 
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t 2 
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w 4 


p t 3 

O 

'\t? 

O 


w 4 


(ID 


where primes denote x derivatives; the four relevant constants are 


Ti = L 2 /D A , T 2 = L 2 /D B) T 3 = l/(kiQ*), T 4 = l/k 2 


and the factor f} resulting from the linearization is 1/(1 + KAo). 

The constants Ti through T4 have the dimension of time, which, even though the problem is steady- 
state, have physical significance since the computed quantities represent transport rates. Unlike the tran- 
sient problem considered in the following section, these parameters always appear as nondimensional 
ratios. 

The first two scales relate to the diffusion of the permeant and carrier species, and the latter two 
relate to the forward and backward kinetic constants. The parameters Da, Dr, ki, and k2 are fixed val- 
ues for the 02-Hb reaction while the membrane length L and the total Hb concentration Q* are subject 
to variation in particular experiments. The second through fourth columns of table 1 show expected val- 
ues of the constants for a representative range of the parameter L using the system parameters from 
Kutchai et al. ( 1970 ). 
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Since the forward chemical time constant T 3 is so fast, the chemical reaction in equation 4 is heavily 
weighed toward the right. In addition, the diffusive time constants are quite large for increasing mem- 
brane thicknesses. Since the reciprocals of these diffusive time constants multiply the space derivatives, 
this term in equation 3 is generally small, which indicates chemical equilibrium. The dual constraint of 
satisfying the equation as well as the boundary conditions causes large deviations from equilibrium at 
the boundaries when the ratio of time constants is large. The details of the solution depend on all four 
time constants, out of which only three independent nondimensional parameters can be formed. The 
Damkohler number, defined as the ratio of the diffusive to chemical time scales, has often been used as a 
representative parameter. The system corresponding to L = 100 p. is characterized as one with a large 
Damkohler number and is referred to as a “fast” reaction. Such linear systems have been thoroughly 
examined in the literature by assuming equilibrium throughout, or by using the methods of asymptotic 
analysis. 

The general solution of equation 1 1 is a sum of four exponentials that determine the spatial dis- 
tribution, each having its own eigenvalue and eigenvector. The particular eigenvalues actually excited 
depends on the boundary conditions. Equation 1 1 has the property that the coefficient matrix possesses 
an eigenvalue equal to zero of rank 2. This simply means that equilibrium solutions (no chemical 
reaction) are admitted with a linear spatial gradient in [O 2 ]. A single parameter that characterizes the 
system is the absolute value of the remaining two eigenvalues which can be calculated from the 
coefficient matrix as 

(n) 

The magnitude of the eigenvalue is computed for a range of membrane thicknesses and is also shown 
in table 1 along with the inverse of the eigenvalue which is a measure of the boundary layer thickness. 
Since one set of eigenvalues is always zero and the other set grows to very large values, this system of 
fast reactions is characterized as a stiff system. There is a large and growing literature concerning stiff 
systems of ordinary differential equations, mostly associated with initial value problems (Seider et al., 
1982). In the case of facilitated transport, the stiffness becomes a problem associated with the boundary 
conditions as well as with the equations themselves, and is manifested by rapid spatial gradients near the 
boundaries. 


SOLUTION OF THE FULL NONLINEAR EQUATIONS 


The nonlinear time-dependent equations are solved by the method of time linearization using an 
implicit numerical algorithm with direct block tridiagonal solvers. In many applications it has been 
found to be advantageous to use direct solvers rather than iterative techniques whenever possible. In the 
unsteady implementation by Ruckenstein, an implicit algorithm was found to be more successful than 
earlier attempts using explicit methods. However, iterative techniques are needed to solve the resulting 
finite difference equations at each time step. Results are presented here using direct solvers for a range 
of cases from slow to very fast chemical reactions, including the temporal evolution of the transport 
process. 
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Most of the developments reported here on the use of implicit numerical methods for stiff systems of 
equations are based on work performed at NASA Ames Research Center in the 1960s on nonequilibrium 
aerodynamics. In Lomax and Bailey (1967) there is a complete discussion of the role of eigenvalues in 
choosing an appropriate numerical scheme for high-speed compressible flow. (It is interesting that the 
local linearization techniques presented in Gonzalez-Femandez (1981) and Moretti (1965) are very simi- 
lar, but are applied in two very different disciplines: biophysics and hypersonic gas dynamics, 
respectively.) 

The method of time linearization involves the simple expedient of expanding the nonlinear term in a 
Taylor series about the solution at a previous time step. The first step is to discretize equation 3 in time 
using Euler forward differences: 


a[o 2 ] to 2 ] - [o 2 ] p 

3t “ At 

(13) 

9[Hb0 2 ] [Hb0 2 ] - [Hb0 2 ] p 

0t = At 

where the subscript p refers to the previous value of the indicated quantity. Next, assemble the resulting 
equation in terms of the vector w defined previously 


w! = F(w, w p ) 

where F is generally a nonlinear function of w at the current and previous time step. The linearized form 
is easily found by expanding F in a first-order vector Taylor series about (jk- ffip)- 


}v' = E(W P , Wp) + (w - Wp)- VF(w p , Wp) 


The net result of this operation is to add two extra terms to the second and fourth rows of the coefficient 
matrix. It also affects the structure of the system by causing the coefficient matrix to vary in both space 
and time. Since no analytical solution is readily available for this linear problem, a numerical approach 
is necessary. In the development of such techniques for multidimensional problems with disparate 
scales, special methods, usually based on stretched coordinates, were developed. These “grid generators” 
are used to cluster the mesh in regions of high parameter variations. A typical application of this 
approach is discussed in Dwyer et al. (1982). In the current problem, with only one space dimension, it 
is useful to consider uniform grids to clearly highlight the resolution required to resolve the boundary 
effects. It is fortunate that the limit solution, that of complete chemical equilibrium, is a known boundary 
point for an infinitely dense grid. 

The final system of equations is: 


y/ = Aw + B 
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where A is the 4x4 coefficient matrix and B is a 4x1 vector. Both A and B depend upon Wp. The matrix 
A is: 


0 




0 


-nw” ! (1 ~ w 3 ) -=r- 
P P T 3 


0 


0 


0 


n M _ 11 

*P t 3 t 4 


0 


O' 


0 


0 0 


The inhomogeneous right-hand side B is: 


0 


T T Tj 

(1 — n)w" — + nw" w 3 rr-wj t- 

*p 1 3 p p I 3 p At 


0 


„ , n l 2 n A 2 A 2 

- (1 -" )w > p tT - nw >p W3 pTT‘ W3 pA7 


The system is solved with second-order central finite differences in vector form. Second-order forward 
and backward formulas are used at the boundaries to avoid including external mesh points. 

After discretization, and incorporation of the boundary conditions, the solution vector w is easily 
computed using the block tridiagonal inversion algorithm described in Isaacson and Keller (1966). As 
mentioned above, one advantage of the matrix formulation is that more reactions and/or species can be 
accommodated by simply increasing the size of the blocks in the tridiagonal system. 


CALCULATED SOLUTIONS 


Membrane thickness is a dominant parameter in facilitated transport. Extremely thick or thin mem- 
branes can be analyzed using equilibrium or linearized methods, respectively. However, many problems 
of practical importance lie in the intermediate range. 

Calculations are shown in figures 2 through 4 for membranes of thickness 1, 10, and 150 microns 
using the system parameters in table 1. Curves are shown for X = [Oil/Co and Y = [Hb02]/Q* from both 
linear and steady-state nonlinear calculations. Also shown are curves of oxyhemoglobin concentration 
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using local equilibrium values of X. At L = 1 p. the distribution of the oxygen permeant is predicted 
quite well using both linear and nonlinear theory. The deviation from complete equilibrium (a linear 
function of X vs. x), is slight, but more apparent is the discrepancy between the equilibrium curve and 
the two computed oxyhemoblobin distributions. There is a significant difference between the linear and 
nonlinear oxyheme distributions near the high pressure side of the membrane. The linear system in 
equation 1 1 is constrained to be symmetrical about x = 0, the midpoint of the membrane, which may 
account for the difference. The discrepancy between the nonlinear Y-curve and the equilibrium curve 
reverses sign at about x = 0.16. This indicates a reaction reversal in equation 6 and subsequent unloading 
of oxygen at the membrane exit. 

Figure 3 shows the same dependent variables, but for a 10-micron-thick membrane. Here the dis- 
crepancy in the computed Y-distributions is obvious and the symmetry constraint seems to drive the 
oxyheme curve to values greater than one — a physically unrealistic result. Another tendency is for the 
oxygen concentration curves (X) to become more distorted, but with the side effect that chemical equi- 
librium is almost complete except for the strong deoxygenation of hemoglobin at the exit. 

Figure 4 shows the most extreme case of L = 150 (i. Here the linear solution is obviously inapplica- 
ble and is not shown. The distribution of oxygen permeant is slightly more distorted, but the most impor- 
tant aspect of this curve is not at all obvious. This is the flux continuity constraint that forces the slope of 
the oxygen curve to match at entrance and exit, since this species is conserved within the membrane. 

The slope (or flux) changes very rapidly within the last one percent of the membrane. The governing 
equations can be manipulated in such a way (Fatt and La Force, 1965) that the flux can be accurately 
computed from the difference in oxyheme concentration across the membrane. The sensitivity of the 
flux (e.g., the difference in Y-values at entrance and exit) to the grid is shown in a sequence of five 
curves with 32, 64, 128, 256, and 512 mesh points. The grid size is shown to have an important effect on 
the computed flux. Figure 5 shows an expanded view of this region and shows clearly how the zero- 
slope boundary condition competes with the oxygen dissociation curve in such a way to induce rapid 
changes near the boundary. The extent of this region is directly related to the eigenvalues from the linear 
analysis as presented in equation 11. 

The tendency toward equilibrium across the membrane is very strong. The only deviation is a very 
thin boundary-layer region where dissociation occurs in order to satisfy the boundary condition. Increas- 
ing mesh density shows an approach to an asymptotic value of Y which is just slightly above the equilib- 
rium value. As mentioned above, the singular perturbation is such that the equilibrium value is always 
available as a lower limit. This effectively fixes the upper limit on the possible flux which can be com- 
puted using algebraic techniques. The size of this boundary- layer region can be approximated from the 
inverse eigenvalue (last column in table 1). 

A calculation showing the unsteady evolution of the permeant is shown in figures 6 to 9 for the case 
L = 10 jx. All calculations are for a 512-point mesh interpolated to 50 points for clarity. Figure 6 depicts 
X as a function of space and time for two different initial concentrations — a linear and an exponential 
distribution. In both cases, chemical equilibrium was assumed in computing the initial values of Y. The 
linear initial distribution in figure 6(a) corresponds to the classical case of passive diffusion through a 
membrane with uniform flux throughout. The distribution of permeant adjusts very rapidly to the imbal- 
ance caused by the boundary condition constraint. This adjustment time is of the order of 0.3 sec which 
corresponds to the slow chemical time scale in table 1. The exponential distribution shown in figure 6(b) 
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is probably a more realistic initial condition since the starting flux vanishes at the membrane exit. Relax- 
ation to steady state takes a little longer. Summary curves of entrance and exit flux are shown in figure 7 
which give a better indication of the approach to steady state. 

These results show that the steady-state solution is independent of initial conditions — an important 
physical confirmation of the mathematical model, since it is by no means assured that long-time 
solutions of nonlinear problems are independent of initial conditions. 

Figures 8(a) and 8(b) are the direct counterparts of figure 6 except for the use of “sigmoid” shaped 
dissociation curves (specifically, curve 2 in figure 1). Aside from a different distribution of species in the 
steady state, the calculated data show a larger transient after the initial start. This transient seems to be 
concentrated in the low-pressure knee of the dissociation curve. The steady distribution of oxygen and 
oxyhemoglobin is shown in figure 9. A comparison with figure 3 shows the following dissociation curve 
effects: 1) Since the sigmoid curve has a high affinity for oxygen, there is very little free hemoglobin at 
the high-pressure boundary; 2) The boundary layer region near the low-pressure end of the membrane is 
wider. This indicates a reduction in stiffness in the governing equations with a less severe requirement 
on the number of grid points. Finally, figure 10 shows the approach to steady state as measured by the 
evolution of the entrance and exit flux. Note the slightly longer settling time and the enhanced excursion 
of the transient when compared with figure 7. 

The most important single parameter in facilitated transport is the overall facilitation factor. Com- 
puted facilitation factors are compared with the experiments of Scholander and the computations of 
Kutchai in figure 1 1 . As mentioned above, the facilitation factor is defined as the ratio of the oxygen 
flux through the membrane to the flux that would exist without a carrier species (Hb in this case) pre- 
sent. The definition used here is one greater than that used by Kutchai et al. (1970). The patterned line in 
figure 1 1 represents the experimental variation of F with variable pC>2* Th e crosses and open circles rep- 
resent computations from Kutchai and the present code, respectively, using the hyperbolic dissociation 
curve. The slight differences are believed to be caused by the differing solubilities used in the two calcu- 
lations. The solid circles represent results using the sigmoidal dissociation curve. Referring to figure 1, 
this curve is a better fit to the dissociation data, as reflected in the better agreement between theory and 
experiment. 

This particular case is a good example of a stiff problem where the facilitation factor is quite close to 
the equilibrium value. This behavior is illustrated in figure 12 for the particular case pC>2 = 25 torr. 
Referring to figure 1, the operating region for the dissociation curve is the range 0 < pC>2 ^ 25. The solid 
curve (dissociation curve 2) shows a more gradual fall in the normalized oxygen concentration at the 
low-pressure boundary when compared to the dashed curve (dissociation curve 1). Figure 13 shows a 
similar comparison for the oxyhemoglobin curves. As expected, the curve obtained using Hill’s expo- 
nent is almost completely saturated at the high-pressure end of the membrane. Finally, figure 14 depicts 
the ratio of the flux of free oxygen to the total (free and oxygenated) flux across the membrane. A value 
of one indicates pure oxygen transport and a value of zero indicates oxyhemoglobin transport only. 
These curves indicate that the physically significant Hill’s curve has a less-severe boundary layer effect 
and that only a small percentage of the oxygen is transported as a free species across the bulk of the 
membrane. This curve also clearly shows the region of active chemical reactions, as indicated by severe 
gradients. The beneficial effect of facilitated transport can be attributed to the fact that a rapid forward 
chemical reaction converts the oxygen to a combined form which is transported across the membrane in 
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this form. At the low-pressure end, the oxygen molecules are “dumped” in a reverse reaction. The strong 
asymmetry caused by the nonlinear interaction seems to be an important aspect of the transport process 
for this stiff problem. 


CONCLUSIONS 


A numerical study of unsteady facilitated transport of oxygen in the nonequilibrium range has been 
extended to model more realistic dissociation processes. A new analysis using implicit numerical 
schemes with direct tridiagonal solvers has been used to solve the governing reaction-diffusion 
equations. 

The most important results from this analysis are: (1) Higher-order reaction models that better match 
oxygen dissociation curve data yield better agreement among theory and experiment for in vitro studies; 
(2) The linearized approximation can give useful information concerning the effect of the system 
parameters, but certain inherent symmetries of the linear theory cause serious errors in the facilitation 
factor as the equations become stiff; (3) The steady-state oxygen flux through a membrane is indepen- 
dent of the initial conditions; (4) The sigmoid oxygen dissociation curve causes almost 100% oxygen 
saturation in the high-pressure side of the membrane and increases the size of the boundary-layer region 
in the low-pressure side. This latter effect relaxes the stringent grid resolution problems using numerical 
methods. 

The computational results reported in this paper were obtained on a VAX minicomputer. The analy- 
sis, however, is quite general and can be extended to greater space dimensions using more powerful 
computer resources to study more complex systems. It is expected that the unsteady aspects will be very 
important for biological systems since the characteristic diffusion time is probably larger than the 
molecular uptake time of a circulating red cell. 
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TABLE 1- Expected Values of Time Constants for Facilitated Transport of Oxygen through 

Hemoglobin Solution 


Da = 1.45x1 0 -5 cm 2 /sec 
Db = 3.00xl0' 7 cm 2 /sec 
ki = 3.18xl0 9 (sec-moles/cc)' 1 
k2 = 40 sec* 1 
Q* = 10x1 O' 6 moles/cc 
p02 =100 torr at entrance 
a = 1.7x1 O' 9 moles/cc/torr 
Co= 1.7x1 O' 7 moles/cc 
Ci = 0.0 


L 

Ti, sec 

T2, sec 

T3, sec 

T4, sec 

X 

\fk 

1 

0.00069 

0.0333 

0.0000314 

0.025 

2.04 

0.4902 

10 

0.0690 

3.33 

0.0000314 

0.025 

20.4 

0.04902 

100 

6.90 

333 

0.0000314 

0.025 

204 

0.00490 

150 

15.5 

750 

0.0000314 

0.025 

306 

0.00327 



OXYGEN PARTIAL PRESSURE p0 2 , torr 


Figure 1.- Dissociation curves for oxygen transport through hemoglobin solution. Curve 1, Second- 
order chemical reaction, p02,o.5 = 7.4 torr; Curve 2, Sigmoid dissociation curve, pO2,0.5 = 7.4 torr; 
Curve 3, Sigmoid dissociation curve for physiological conditions, pO2,0.5 = 26 torr. Symbols, 
experiment (Hemmingsen, 1960). 
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DISTANCE ACROSS MEMBRANE, x 


Figure 2 - Distribution of X = [02]/Co and Y = [HbC>2]/Q* across the membrane. Second-order 

chemical reaction, L = 1 micron, p02 = 100 torr at entrance, 0 at exit. , nonlinear calculation; 

, linear calculation; equilibrium curve Y(X). 
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Figure 3 - Distribution of X = [C> 2 ]/Co and Y = [HbC>2]/Q* across the membrane. Second-order chemi- 
cal reaction, L = 10 microns, pC >2 = 100 tore at entrance, 0 at exit. , nonlinear calculation; — 

, linear calculation; equilibrium curve Y(X). 
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DISTANCE ACROSS MEMBRANE, x 

Figure 4 - Distribution of X = [O 2 VC 0 and Y = [HbC> 2 ]/Q* across the membrane. Second-order chemi- 
cal reaction, L = 150 microns, pC>2 = 100 torr at entrance, 0 at exit. , nonlinear calculation; - 

, equilibrium curve Y(X). Curves of Y shown for 5 different computational grids. 
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Figure 5 - Same conditions as figure 4 showing low-pressure side of membrane. 
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Figure 6 - Space-time distribution of oxygen concentration. Second-order chemical reaction, 

L = 10 microns, p02 = 100 torr at entrance, 0 at exit, a) Constant flux (linear) initial distribution, 
b) variable flux (exponential) initial distribution. 
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Figure 7 - Transient flux (normalized by flux without chemical reaction, Fo) at entrance (ENT) and exit 

(EX) of membrane. Same parameters as figure 6. , linear initial distribution; , 

exponential initial distribution. 


22 




- [O2I /Cq 


1.0 



Figure 8 - Space-time distribution of oxygen concentration. Sigmoid dissociation curve, L = 10 microns, 
pC >2 =100 torr at entrance, 0 at exit, a) Constant flux (linear) initial distribution, b) variable flux 
(exponential) initial distribution. 
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Figure 9 - Distribution of X = [C> 2 ]/Co and Y = [Hb02]/Q* across the membrane. Sigmoid dissociation 

curve, L = 10 microns, pC >2 = 100 torr at entrance, 0 at exit. , nonlinear calculation; , 

linear calculation; equilibrium curve Y(X). 
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Figure 10 - Transient flux (normalized by flux without chemical reaction, Fo) at entrance (ENT) and exit 

(EX) of membrane. Same parameters as figure 8. , linear initial distribution; , 

exponential initial distribution. 
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Figure 11- Variation of steady-state facilitation factor with oxygen partial pressure at membrane 
entrance. L = 150 microns, other parameters as in table 1. = = = = = =, experiment (Scholander); 
solid symbols, sigmoid dissociation curve; open symbols, second-order dissociation curve; crosses, 
computation (Kutchai) using second-order dissociation curve. 
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Figure 12 - Comparison of X = [C>2]/Co distributions across a membrane. pQz = 25 torr at entrance. 
, sigmoid dissociation curve; , second-order dissociation curve. 
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Figure 13 — Comparison of Y = [HbC> 2 ]/Q* distributions across a membrane. pC >2 = 25 torr at entrance. 
, sigmoid dissociation curve; , second-order dissociation curve. 
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Figure 14.— Transport of free oxygen fraction across a membrane. pC>2 = 25 torr at entrance. , 

sigmoid dissociation curve; , second-order dissociation curve. 
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